      Program ex3
******************************************************************************
*                                                                            *
*  This program is to illustrate the use of 2D and Multiple 1D filters       *
*  We use here DFIR2D and DFIRM1D modules from LIBCONV                       *
*                                                                            *
******************************************************************************
*                                                                            *
*  The 1D filter in the X direction will be (1/2,-1.,1/2)                    *
*                               --- Approximation of second derivative       *
*  In the Y direction the 1D filter will be (-1/2,0.,1/2)                    *
*                               --- Approximation of first  derivative       *
*  The equivalent 2D filter will be:                                         *
*                               -1/4   0   +1/4                              *
*                               +1/2   0   -1/2                              *
*                               -1/4   0   +1/4                              *
*                               ------ Y ----->                              *
*  This is the approximation of the triple derivative D()/Dxxy               *
*                                                                            *
******************************************************************************
      parameter (HX=5, HY=5, NX=(2*HX+1), NY=(2*HY+1))
      parameter (PI = 3.14159265358979323846)
      double precision input(-HX:HX,-HY:HY), fil2d(-1:1,-1:1)
      double precision filx(-1:1), fily(-1:1), temp(-HX:HX,-HY:HY)
      double precision out1(-HX:HX,-HY:HY), out2(-HX:HX,-HY:HY)
      double precision scale_x, scale_y
      data fil2d/-.25,0.5,-.25,0.0,0.0,0.0,.25,-.5,.25/
      data filx /0.5,-1.,0.5/,        fily/-.5,0.0,0.5/
* Init Input
      scale_x = PI/dfloat(HX)
      scale_y = PI/dfloat(HY)
      do j = -HY,HY
        do i = -HX,HX
          input(i,j) = cos(i*scale_x)*(1.-dfloat(ABS(i))/dfloat(HX))
     $               * cos(j*scale_y)*(1.-dfloat(ABS(j))/dfloat(HY))
        end do
      end do

*  2D FIR filter:
      call dfir2d( input,1,NX,-HX,NX,-HY,NY,
     $             fil2d,1,3, -1,  3, -1, 3,
     $             out1, 1,NX,-HX,NX,-HY,NY, 1.0d0, 0.0d0)

*  Multiple 1D FIR filter in the X direction:
      call dfirm1d( input,1,NX,-HX,NX, NY,
     $              filx, 1,   -1, 3,
     $              temp, 1,NX,-HX,NX,       1.0d0, 0.0d0)

*  Multiple 1D FIR filter in the Y direction:
      call dfirm1d( temp,NX,1,-HY,NY, NX,
     $              fily, 1,   -1, 3,
     $              out2,NX,1,-HY,NY,        1.0d0, 0.0d0)

*  Print out the input and the two results
      write(6,11)
      write(6,12) ((input(i,j), j=-HY,HY), i=-HX,HX)
      write(6,13)
      write(6,12) ((out1(i,j), j=-HY,HY), i=-HX,HX)
      write(6,14)
      write(6,12) ((temp(i,j), j=-HY,HY), i=-HX,HX)
      write(6,15)
      write(6,12) ((out2(i,j), j=-HY,HY), i=-HX,HX)
 11   format('Input : ')
 12   format(11f7.3)
 13   format(/'After application of 2D FIR filter:')
 14   format(/'After multiple 1D FIR filters in the X Direction:')
 15   format(/'After final multiple 1D FIR filters in the Y Direction:')
      stop
      end
